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ABSTRACT 

We have modeled the emission from dust in pre-protostellar cores, including a self-consistent calculation 
of the temperature distribution for each input density distribution. Model density distributions include 
Bonnor-Ebert spheres and power laws. The Bonnor-Ebert spheres fit the data well for all three cores 
we have modeled. The dust temperatures decline to very low values (T^ ~ 7 K) in the centers of 
these cores, strongly affecting the dust emission. Compared to earlier models that assume constant dust 
temperatures, our models indicate higher central densities and smaller regions of relatively constant 
density. Indeed, for L1544, a power-law density distribution, similar to that of a singular, isothermal 
sphere, cannot be ruled out. For the three sources modeled herein, there seems to be a sequence 
of increasing central condensation, from L1512 to L1689B to L1544. The two denser cores, L1689B 
and L1544, have spectroscopic evidence for contraction, suggesting an evolutionary sequence for pre- 
protostellar cores. 



Subject headings: stars: formation — ISM: dust, extinction — ISM: clouds 
(L1689B, L1544, L1512) 



ISM: individual 



1. INTRODUCTION 

Molecular cloud cores that have not yet formed a star 
provide the best opportunity to determine the initial con- 
ditions for star formation. The density distribution as a 
function of radius in such starless cores is a strong dis- 
criminator between theoretical models, and it provides an 
essential tool for understanding the process of collapse and 
star formation that follows (Andre, Ward-Thompson, & 
Barsony 2000). However, the earliest evolution of a cloud 
core is very poorly constrained, and its evolutionary track 
is unclear. 

For those more evolved objects that possess infrared 
sources, the sequence is well described by the empirical 
evolutionary sequence of Lada (1987), which is supported 
by theoretical modeling of the spectral energy distribution 
(SED) (Adams, Lada, & Shu 1987). In this scheme, the 
well-known classification (i.e. Class I, II and III) is based 
on the near-infrared spectral index and is interpreted in 
terms of the sources becoming progressively less embed- 
ded in circumstellar dust. The scheme was extended to 
younger, more deeply embedded objects with the Class 

classification of Andre, Ward-Thompson, & Barsony 
(1993). These objects possess a highly enshrouded IRAS 
source (so that their SEDs peak longward of 100 /im); ac- 
cretion models suggest that the mass of the central object 
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is less than or similar to that of the collapsing envelope. 
Once a protostellar core has formed, the density distribu- 
tions are strongly centrally peaked and the evolution of 
the objects has a "singular" nature in that the dynam- 
ics are primarily determined by the characteristics of the 
innermost regions of the cores. 

At earlier stages of evolution, starless cores, defined as 
dense cores of gas that do not possess an IRAS source (My- 
ers & Benson 1983; Benson & Myers 1989), are potential 
sites of star-formation. An important breakthrough in the 
study of these objects was the detection of dust continuum 
emission at millimeter and submillimeter wavelengths by 
Ward-Thompson et al. (1994), who dubbed the subset of 
starless cores that possess such emission "pre-protostellar 
cores" (hereafter PPCs). 

However, unlike the Class O-III sources described above, 
the classification of a source as a PPC can have a much 
broader interpretation: the density distribution in these 
objects could be anything between uniform and strongly 
centrally peaked. They may or may not be gravitation- 
ally bound. Moreover, the dynamical and evolutionary 
status of PPCs is poorly defined; they may be quasi- 
statically contracting (perhaps supported by non-thermal 
pressures), or in an early stage of dynamic collapse. Al- 
ternatively, it is possible that a combination of magnetic, 
turbulent and other effects may inhibit collapse or even 
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suppress it completely. 

It is crucially important to establish the physical con- 
ditions within the extended envelopes of PPCs since they 
will determine most of the early evolutionary history of 
the sources as they evolve towards main-sequence stars. 
In particular, the initial conditions are highly important in 
defining the collapse dynamics, the likely mass of the star 
(and globally, the initial mass function), the timescales 
for evolution, and the likelihood of detection of sources in 
different evolutionary epochs. 

Difi'erent initial conditions lead to very difi'erent collapse 
dynamics. Collapse from a singular isothermal sphere 
leads to the "inside-out" collapse model with a constant 
mass accretion rate (Shu 1977). If collapse begins before 
the density distribution relaxes fully to a singular isother- 
mal sphere, there will be an initial period of much higher 
mass accretion rate, followed by a declining rate of accre- 
tion (Foster & Chevalier 1993; Hcnriksen, Andre, & Bon- 
temps 1997). Such an epoch has been identified with 
the Class objects (Andre, Ward-Thompson, & Barsony 
2000). Collapse from a logatropic sphere, which produces 
mass accretion rates that increase with time, has also been 
suggested (McLaughlin & Pudritz 1997). To address these 
issues, it is necessary to understand the density and tem- 
perature distributions within the cores and to correlate 
them to the kinematics. In this paper we address the first 
of these issues. 

In the original studies of PPCs (Ward-Thompson et al. 
1994; Ward-Thompson, Motte, & Andre 1999), the ra- 
dial distribution of the continuum emission was found to 
be relatively flat toward the center of the PPCs; this in- 
tensity distribution was interpreted purely in terms of a 
density distribution that is relatively flat within the in- 
ner core and falls off as an approximate power law at 
larger separations. These authors, and others, have in- 
voked such density distributions as proof of the significance 
of non-thermal (magnetic and/or turbulent) pressure sup- 
port in these objects. However, Ward-Thompson, Motte, 
& Andre (1999) pointed out that the source statistics 
and, where available, the observed magnetic field strengths 
suggest a quantitative disagreement with the evolution- 
ary timescales in the published ambipolar diffusion models 
for the quasi-static contraction of magnetized cloud cores. 
This conclusion has been questioned by Ciolek & Basu 
(2001). 

These conclusions arc important in constraining the col- 
lapse process, but previous interpretations of the data as- 
sumed that the dust temperature is uniform throughout 
the cores. We will argue that the interpretations of the 
observational data are profoundly affected by this assump- 
tion. Early calculations for dust clouds heated by the ex- 
ternal radiation field found that the temperature declined 
to very low values in the interior (Leung 1975); more re- 
cently Zucconi et al. (2001) have applied a semi-analytic 
technique and come to similar conclusions. 

In this paper, we constrain the density and temperature 
profiles of a small sample of PPCs, utilizing high quality 
observational data and a dust continuum radiative transfer 
code (Egan, Lcimg, & Spagna 1988). We consider a range 
of physical models for the density distribution, calculating 
the dust temperature distribution in each ease in a self- 
consistent way. The model clouds are then "observed" 



to mimic the observational technique (including the ef- 
fects of the telescope beam power pattern and chopping) 
in order to compare to actual observations. Although sev- 
eral sources of observational data are considered, the pri- 
mary observational comparison is to our recent 450/im and 
850/xm data obtained with the SCUBA on the James Clerk 
MaxweU Telescope (JCMT) (Shirley et al. 2000). 

2. PHYSICAL MODELS 

We consider the simplest physical situation of thermally 
supported, gravitationally bound, spherically symmetric 
clouds. While all three cores modeled here have some de- 
gree of turbulence, it is subsonic in L1512 and L1544 or 
barely sonic in L1689B; the effects of turbulence should not 
be dramatic in these very quiescent cores. The solution for 
the thermally-supported, isothermal cloud is the Bonnor- 
Ebert sphere, which can be generalized to account for a 
temperature gradient (see below). In addition to the ther- 
mally supported Bonnor-Ebert spheres, we have also con- 
sidered a singular sphere, for which p oc r~^, normalized 
to match a Bonnor-Ebert model at the outer radius. Mag- 
netically supported cores have also been suggested (Ciolek 
& Mouschovias 1994; Ciolek & Basu 2000). The latter 
authors have presented density profiles that are specifi- 
cally appropriate to L1544. Although a comparison with 
these models would be of great interest, the inclusion of 
magnetic pressures necessarily implies significant devia- 
tions from spherical symmetry - in fact, highly flattened 
disk-like structures are predicted. After some attempts to 
capture the essence of these models in the one dimensional 
models we consider here, it became clear that a fully two 
dimensional treatment is required. We defer discussion 
of these magnetic models to future work but discuss the 
possible effects of asphericity on our modeling in §8. 

We have considered two types of hydrostatic pressure- 
balanced cloud: 

1. "Classical" isothermal Bonnor-Ebert spheres, in 

which (dT/dr) = 0; and 

2. Modifled Bonner-Ebert spheres with a gradient in 
the kinetic temperature. 

If the gravitationally bound clouds are non-magnetic, 
arc not subject to any large-scale systemic velocity flows 
(rotation, infall, or outflow), and are purely supported by 
thermal pressure, then the equations of hydrostatic equi- 
librium are as follows: 

dP{r) GM{r)p{r) 



dr 

dM{r) 

dr 



= 47rr^p(r), 



(1) 
(2) 



where r is the radius, P{r) and p{r) are the pressure and 
density respectively at r, and M(r) is the total mass en- 
closed within r. If the clouds are isothermal, the equation 
of state is 

P = a'^p, (3) 

where a is the isothermal sound speed (a^ = kT / pmn, k is 
the Boltzmann constant, /i is the average molecular mass 
in the gas, and tuh is the mass of the hydrogen atom). 

With an appropriate change of variables to a dimen- 
sionless form, D = p/ pc and ^ = [r / a)-/^KUpl, where Pc 
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is the central density, these equations yield a single solu- 
tion, with the maximum value of ^ (^max.) as the only free 
parameter (cf. Foster & Chevalier, 1993), which is deter- 
mined by pressure balance with the surrounding medium. 
If £,max. > 6.451, which corresponds to a eenter-to-edge 
density contrast of po/pc > 14.3, then the sphere is un- 
stable to perturbations and subsequent gravitational col- 
lapse. The marginal case is the "critical Bonnor-Ebert 
sphere" (Ebert 1955, Bonnor 1956). The solution corre- 
sponds to a family of spheres of differing central density 
[cf. Fig. 1 of Shu (1977)]. The density structure of the 
outer envelopes is closely approximated by p oc r~^, but 
the density approaches a constant value at small r. For a 
given temperature, higher values of pc yield models with 
smaller cores. The most extrcinic; form of the; unstable 
Bonnor-Ebert sphere has an infinite central density (cor- 
responding to ^max. = oo)- This is the singular isothermal 
sphere, in which p oc throughout. 

The modified Bonnor-Ebert spheres include the effects 
of radial temperature variations. In these cases, 



dP 

dr 



= k 



dT 
dr 



T- 



^dn 
dr 



(4) 



and, using p = firnHn, where n is the number density, 
equation 2 becomes 

dM{r) 



dr 



{A'Kp,mH)nr 



and, from hydrostatic equilibrium, 

dn / p,mHG\ M{r)n fn\dT 
dr 



Tr^ 



T) dr 



(5) 
(6) 



In all of our models the mean molecular mass is /x = 2.3, 
consistent with n = n(H2) + n(He). The inner radius 
(ri) is not constrained by observations and is simply set 
small enough (25 AU) that it does not affect the results, as 
shown by tests that used = 50 AU. The free parameters 
in the isothermal Bonnor-Ebert models are thus the cen- 
tral density (ric), the outer radius (r,,) and the gas kinetic 
temperature {Tk)- Most models have an outer radius of 
3 X 10* AU, but some have an Vo of 1.5 x 10* or 6 x 10* 
AU. 

For most of our models (see Table 1) we have used the 
isothermal Bonnor-Ebert configurations with Tk = 10 K. 
The density distribution from a sequence of these mod- 
els (with different central densities, ric) is shown in Figure 
la. A few isothermal models were computed for differ- 
ent values of Tk to examine the sensitivity of the density 
structure to the temperature (Fig. lb). The main effect 
of changing the assumed Tk is to change the size of the 
central core and the mass (Table 1). 

For all cases with tt-c > 3 x 10^ cm~^, the density con- 
trast between the edge and the center of the cores exceeds 
14.3, implying that they are unstable to collapse. For the 
isothermal Bonnor-Ebert spheres of fixed temperature and 
outer radius, the enclosed mass is only slightly sensitive to 
the central density (nc). Since the photometric fluxes are 
largely determined by the masses within the beams, ob- 
servations of the spectral energy distribution (SED) alone 
cannot place very strong constraints on the central den- 
sities of Bonnor-Ebert spheres. However, by modeling 
the radial intensity distribution it is possible to disc-rimi- 
nate between models with different central densities, which 
have different density distributions. 



In order to generate the non-isothermal Bonnor-Ebert 
spheres it is necessary to pre-define the temperature pro- 
file. In practice, this is achieved by using an appropri- 
ately defined isothermal Bonnor-Ebert model as the input 
for the radiative transfer code to generate the dust tem- 
perature profile (see §3). This Td{r) is then used to re- 
calculate the density profiles in the non-isothermal model. 
Obviously, the corrected density profile could then be used 
to generate a second iteration for the temperature profile 
etc., but the magnitude of the corrections in the second 
and higher order iterations is sufficiently small for them to 
be neglected. The correction for non-isothermality results 
in a smaller core and an initially steeper density profile 
(Fig. Ic), caused by the second term in equation 6, with 
dT/dr > 0. Beyond radii of about 5000 AU, the first term 
in equation 6 dominates, and the smaller enclosed mass 
(caused by the smaller core) leads to a shallower density 
profile. The total mass is larger because the mass is dom- 
inated by the outer regions. The shallower density profile 
in the outer envelope decreases the center to edge density 
contrast (e.g. in the case of ric = 10^ cm~^, = 3 x 10* 
AU, by a factor of ~1.8). 

Note that the radiative transfer code computes the ra- 
dial profile of the dust temperature, Td{r), whereas the 
equilibrium configuration is determined by the gas kinetic 
temperature, TK{r). While efficient gas-dust coupling 
forces Tk to equal Td at high densities, at densities be- 
low about 1 X 10* to 3 X 10* cm^^ Tk ^ Td (Takahashi, 
Silk, & Hollenbach 1983; Doty & Neufeld 1997). Al- 
though these models are a theoretical improvement on the 
isothermal Bonnor-Ebert sphere, a full calculation of the 
gas energetics, including dust coupling, would be needed 
to make fully self-consistent models. 

All the models of the density profiles are tabulated in 
Table 1. All densities are specified as the number density 
of particles. The notation tBE implies a non- isothermal 
Bonnor-Ebert model in which the density structure was 
computed assuming that TK{r) = Td{r), where Tdir) was 
computed from the initial isothermal BE model with the 
same ric (model number denoted by the superscript). We 
again emphasize the distinction between Tk and Td by us- 
ing the term isothermal only to refer to the assumption 
about TK{r) used to compute the density distribution. 
The dust temperature in these isothermal BE models is 
not constant (§3). 

3. INTERSTELLAR RADIATION FIELD 

The radiative transport code of Egan, Leung, & Spagna 
(1988), modified to use an arbitrary density distribution, 
computed Td(r) for each physical model. We first explored 
the dependence of the temperature profiles on the assumed 
form of the interstellar radiation field (ISRF). Figure 2 
shows a comparison of two estimates of the ISRF. Pre- 
vious calculations with the radiative transfer code (e.g., 
Leung (1975); Zhou et al. (1990)) have used an ISRF 
similar to that of Mathis, Mezger, & Panagia (1983), 
supplemented by a blackbody for the cosmic background 
radiation (hereafter "MMP"). More recent analyses, us- 
ing the COBE data, indicate a stronger ISRF in the in- 
frared (Black 1994). Significant departures from the 
"MMP" ISRF occur between 5 and 400 p.m. We have 
modified the Black (1994) ISRF at the shorter ultravi- 
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olet wavelengths (A < 0.36 fim) using equations given 
by van Dishoeck (1988), which reproduce the ISRF of 
Drains (1978). The differences between the "MMP" and 
the Black-Draine models of the ISRF are quite substantial, 
reaching a factor of 1.8 in parts of the ultraviolet, and a 
factor of 13 in parts of the infrared. 

Figure 3 compares temperature profiles using the two 
different models of the ISRF, with the same physical 
model: a Bonnor-Ebert sphere with a central density of 
1 X 10® cm~'^ and uniform kinetic temperature, Tk = 10 
K. The differences in Td{r) are small (< 1.2 K), but the 
biggest differences are at the cold center of the cloud, 
where a small change in can potentially have a major 
impact on the emission. Comparison of the emitted fluxes 
calculated for different wavelengths and beams indicate 
modest changes (less than 15%) for A > 170 /xm, but big- 
ger changes at shorter wavelengths. In our modeling, we 
use the Black-Draine model as the standard representation 
of the ISRF. 

4. DUST OPACITIES 

Dust opacities in molecular clouds clearly differ from 
those in the general interstellar medium. Observations of 
regions forming massive stars (van der Tak et al. 1999, 
2000) have been well-matched by a set of opacities cal- 
culated for grains that have grown by coagulation and 
accretion of thin ice mantles for 10^ years at a density 
of 10^ cm~^, listed in column 5 of Ossenkopf & Henning 
(1994) (hereafter 0H5). This model of the dust opacity is 
therefore probably appropriate to the cold, quiescent cloud 
cores that are the subject of this study. We adopt these 
opacities as the standard model and explore the effect of 
using different opacities. As discussed below, the excellent 
agreement between the observational data and our physi- 
cal models strongly supports the validity of these opacity 
laws. 

As a specific alternative to 0II5, we consider the opac- 
ities in column 2 of Ossenkopf & Henning (1994) (here- 
after 0H2), which include coagulation but lack ice man- 
tles. These may be more appropriate in regions where 
young stars have heated the grains enough to remove the 
mantles. The model opacities are larger at 1300 yum than 
in the 0II5 grain models by a factor of about 2.2 (Fig. 2). 
Both the 0H5 and the 0H2 dust models result in consid- 
erably higher opacities at long wavelengths (by a factor 
of 6.4 for 0II2 at 1300 /zm) than for models in which the 
grains have not undergone coagulation (e.g., column 1 of 
Ossenkopf & Henning (1994). 

However, despite these differences, we find that the ef- 
fect of the opacity law on Td{r) is small; the largest differ- 
ence is at the inner radius, where 0H2 opacities give a 
that is lower by 0.5 K than that obtained with the 0H5 
opacity model. 

5. ALTERNATIVE HEATING MECHANISMS 

Because the very low temperatures of the dust grains 
in the centers of dense cores have a major impact on the 
interpretation of the submillimeter emission (§6 to §7), it 
is important to consider all other possible heating mecha- 
nisms in addition to the absorption of the ambient ISRF. 
The primary source of energy deep inside opaque cloud 
cores is provided by cosmic rays. We consider direct depo- 



sition of energy in dust grains by cosmic rays, absorption of 
ultraviolet photons produced by secondary cosmic ray ex- 
citation of H2 molecules, and, finally, energy transfer from 
possibly warmer gas heated by the cosmic rays. For the 
purpose of a comparative study, we consider a representa- 
tive grain of radius 0.1 /im, absorption properties given by 
0H5 dust, and a cosmic ray ionization rate ^ = 1 x 10"""^^ 
s~^. This value is near the maximum found in recent inves- 
tigations (de Boisanger, Helmich, & van Dishoeck 1996; 
Caselli ct al. 1998). The most recent estimates indicate 
C = (2.6±1.8) X 10^^^ s-\ based on modeling of H13C0+ 
observations (van der Tak & van Dishoeck 2000), which 
is in good agreement with direct measurements by the 
Voyager and Pioneer spacecraft (Webber 1998). Con- 
sequently, the value that wc have adopted should provide 
a strong upper limit to the cosmic ray heating. In each 
case, we compare the heat input to a grain to the radia- 
tive cooling by a grain at a temperature of 5 K, somewhat 
lower than the value of reached at the center of the 
cores. Because the dust in the code is in radiative equilib- 
rium, the cooling rate sets a lower limit on the heating rate 
due to the ISRF. At this temperature the power emitted 
by our representative grain, obtained by integrating over 
the product of the Planck function and the emission cross 
section for 0H5 dust, is 4.0 x 10"^^ erg s"^ 

Firstly, we consider the direct cosmic ray heating of dust 
particles. As discussed by Greenberg (1991), the energy 
input is maximized if we assume that all cosmic ray ion- 
izations are caused by 1 MeV protons. Then, our value of 
C implies a flux of protons of 4.8 cm~^s~^, each of which 
can deposit 6.13 keV in the representative grain described 
above (Greenberg 1991). This results in an energy de- 
position rate per grain of 2.1 x 10"^'' erg s~"^, 190 times 
smaller than the radiative heating. 

Secondly, we consider the effect of ultraviolet photons 
created following the cosmic ray ionization of H2. This 
process yields energetic electrons that are capable of col- 
lisionally exciting electronic states of H2; which decay 
producing a spectrum of ultraviolet photons (Prasad & 
Tarafdar 1983). Estimates of the flux of these photons 
range from 1.4 x 10^ (Prasad & Tarafdar 1983) to 1 x lO"* 
cm~^s^^ (Greenberg & Li 1996). A calculation of partic- 
ular relevance is that of Cecchi-Pestellini & Aiello (1992), 
who include the effects of a density gradient and internal 
extinction of the photons. They consider a centrally con- 
densed cloud with properties similar to those that we infer 
for pre-protostellar cores and find a maximum photon fiux 
near the center of of 'I'^ = 1 x 10^ photons cm~^ s~^, for 
^ = 4 X 10~^^ s~^. Scaling to our adopted value of ( gives 
$1/ = 2.5 X lO''' photons cm~^ s~^ as a maximum value. 
If all of these photons are at a wavelength of 120 nm [this 
probably overestimates the average photon energy, (Gre- 
del et al. 1989)], then the energy deposition per grain is 

^^E^QaTra^ < 2.5 x 10^x1.7 x 10-"x5.1 x 10"^° = 2.2 x 10 

(7) 

which is 18 times smaller than the heating rate due to the 
external ISRF. 

Finally, we consider heating of dust grains by collisions 
with warmer gas particles. The energy deposition rate de- 
pends on the gas temperature, which is poorly constrained. 
However, we can put a limit to the total heating caused by 
cosmic ray heating of the gas by assuming that all of the 
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energy input to the gas is transferred to the dust. This is 
clearly a gross upper limit since much of the energy will 
be radiated away in molecular emission. The volume rate 
of heating by cosmic rays (Goldsmith & Langer 1978) is 

TcB. = 3.2 X 10-27 (_-i_^)n(H2)erg cm-^ s-\ (8) 
Assuming that 

= n(H2)/3.9 X 10", (9) 

where n(H2) and na are the number densities of H2 
molecules and dust grains, respectively, the upper limit 
of the heating rate per grain is 1.3 x 10^^"'' erg s^^; still a 
factor of 3 less than the heating by the ISRF. In princi- 
ple, we could also consider the heating of the gas by the 
ISRF in the outer parts of the core, but our main concern 
is whether or not the central regions of the cores arc as 
cold as our models predict. In any case, the density in the 
outer regions is too low for the gas and dust temperatures 
to be well coupled. We conclude that heating by the ISRF 
dominates other possible heat sources for dust grains, even 
deep inside an opaque core. 

Recent calculations of gas and dust energetics, includ- 
ing the effects of depletion on the coolants, support the 
conclusion that the gas cannot substantially raise the dust 
temperature (Goldsmith 2001). 

6. MODELING OBSERVATIONS 

The primary observational data that the models must 
reproduce are the radial profiles of intensity at 450 and 
850 /Ltm, as presented by Shirley et al. (2000). For some 
sources, intensity profiles at 1300 /zm obtained with the 
IRAM 30-m telescope also exist. We model the radial 
intensity distribution, /""'''"(t) = I^(^b)/I^{bo), as a func- 
tion of impact parameter (b), normalized to the innermost 
impact parameter, corresponding to one quarter of the 
beamwidth. This measure is thus sensitive to the nor- 
malized density and temperature distributions, while the 
SED is sensitive to the ISRF and a quantity that is propor- 
tional to the product of total mass within the beam and 
the opacity. Thus, modeling the normalized radial inten- 
sity distribution and the SED separately provide roughly 
orthogonal constraints on different model variables. Pho- 
tometric data, summarized by Shirley et al. (2000), can 
constrain the ISRF, mass, and model opacities. Of course, 
photometric data obtained with different beams at the 
same wavelength retain some sensitivity to the density and 
temperature distributions. 

To model correctly the radial intensity distribution, one 
must calculate the emission from a model cloud, with a 
self-consistent Td(r), convolve the emission with the ob- 
served beam shape, and "chop" the model in the same way 
as the data were chopped during the observations. The cal- 
culation of Td{r) is done with the code of Egan, Leung, & 
Spagna (1988). Figure 4 shows that more centrally con- 
densed cores (higher nc) result in lower dust temperatures 
in the interior (as expected), but somewhat higher dust 
temperatures at larger radii until all models converge to 
about the same value at the edge of the cloud. 

A second code then uses these density and temperature 
profiles [n(r) and 7d(r)] to calculate the angular depen- 
dence of the emission at specific wavelengths, convolves 
with the observed beam, and uses a numerical simulation 



of the chopping that was used in obtaining the data to 
produce a predicted intensity profile, normalized in the 
same way as the data (-?^o™(^))- The simulation of chop- 
ping in a one dimensional code cannot replicate exactly 
the observations at large impact parameters. We limit our 
comparison of models to data to ranges of b where the pre- 
dictions are insensitive to the exact model of chopping. In 
addition, the ratio of the normalized intensities at the two 
wavelengths [12^™ {b) / I^^q"^ (b)] is computed to examine 
possible variations in the spectral index with radius. Fig- 
ure 4 shows how the predicted intensity distributions are 
affected by the assumed form of the density distribution. 
Although the more condensed models lead to more rapid 
declines of the density with radius, beam convolution and 
chopping, as shown in Figure 4c; and 4d, makes the differ- 
ence less striking than one might expect [see also Fig. 1 of 
Shirley et al. (2000)]. 

For the 1300 ^m data obtained with the 30-m IRAM 
telescope, the beam power pattern is available from the 
literature, but the data cannot be completely modeled be- 
cause they were obtained with multiple chopper throws 
and restored. To bracket the range, we ran models with 
no chopping and models with a 120" chop, based on an 
estimate of where the IRAM data begin to lose sensitiv- 
ity (P. Andre, pers. comm.). Clearly, some loss of signal 
at large impact parameters has occurred because the im- 
chopped models do not reproduce the declines seen in the 
multiply-chopped data, but chopping at 120" may overes- 
timate the effect. 

Modeling the SED allows tests of different aspects of the 
models: the ISRF, the assumed dust opacities, and the 
mass within the beams. In particular, far-infrared data 
constrain the ISRF, while submillimeter data constrain 
the product of mass and submillimeter opacity. Photo- 
metric data are quite variable in quality; most references 
do not provide a beam shape and the data often repre- 
sent the emission integrated over a map. In these cases, 
we assume a Gaussian beam. In addition, calibration un- 
certainties are hard to quantify. Far-infrared data from 
ISO exist (Ward-Thompson, Andre, & Motte 1998) and 
(Ward-Thompson & Andre 1999), but the flux density es- 
timates were being revised to correct for background sub- 
traction. While this paper was being refereed, the new 
measurements became available (Ward- Thompson et al. 
2001). The new results require an ISRF lower than used 
in our models, but the conclusions about the best density 
model are unaffected. Consequently, only a few new mod- 
els were calculated for each source to illustrate the effects 
of the lower ISRF. These new results strongly support our 
conclusions that the dust temperatures are very low in the 
centers of these cores 

The agreement between the model and the data can be 
quantified in terms of the reduced chi-squared 

Xl = {I~{b) - IZTib))MyN (10) 

where a is the uncertainty in the data and N is th(^ mim- 
ber of data points; only points spaced by a full beam are 
used in computing Xr to avoid introducing correlations. 
Because the Xr measure gives great weight to a few dis- 
crepant points and/or points with small uncertainties, we 
have also computed the mean absolute deviation, as used 
in robust estimation: 

(<5) = |/— (6)-/— (6)|/Ar. (11) 
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The statistics were computed separately for each obser- 
vation. Generally, the computations were carried out to 
angular distances of about 60", beyond which inability to 
simulate chopping exactly make the comparison dubious 
for the SCUBA data. A similar calculation with the flux 
densities (S^) was used for comparing the observed and 
model SEDs. In Tables 2 to 4, we report the sum of the 
xl for the 450 and 850 /um data, as well as {6), but only 
Xr for the comparison to 1300 /xm data and the SED. The 
sum of Xr for the SCUBA data is always dominated by the 
850 /Ltm data because of the smaller uncertainties, whereas 
the (i5) weights the 450 /im and 850 /zm data more equally. 
For the 1300 /xm data, two values for Xr given in Ta- 
bles 2-4, the chopped value first, and the unchopped value 
second. In the figures, the chopped prediction for 1300 /im 
is shown as a dashed line, whereas the solid line indicates 
unchopped models. 

The x^ values for the SED in most tabic entries do not 
include the new far-infrared data of Ward-Thompson et 
al. (2001); two models for each source, at the bottom of 
each table, do include these data in the Xr calculation. 

7. MODELLING INDIVIDUAL SOURCES 

7.1. L1689B 

For L1689B, we assume a distance of 125 pc, rather 
than the traditional 160 pc, based on improved determi- 
nations by de Gcus. dc Zccuw, & Lub (1989), and Kmidc 
& Hog (1998), which arc supported by the distance to 
upper Scorpius (de Zeeuw ct al. 1999). The SCUBA 
data on this source (Shirley et al. 2000) have fairly good 
signal-to-noise and the source shape as projected onto the 
plane of the sky is reasonably circular, away from the 
central region. The higher contour levels are somewhat 
elongated and may even have weak multiple peaks. Thus, 
our spherical models can only be approximations to the 
actual source. In addition to the SCUBA data, we con- 
sider the IRAM 1300 fim map (Andre, Ward-Thompson, 
& Motte 1996). As stated above, because these data 
were taken with multiple chopper throws and restored, we 
cannot simulate the effects of chopping as we can for the 
SCUBA data; thus we weight the agreement with the 1300 
/im data much less in deciding which models are best. 

For most of the models of L1689B, whose results are 
summarized in Table 2, we adopt the simple, isothermal 
Bonnor-Ebert description of the density distribution. In 
these models, the primary variable is the degree of cen- 
tral condensation, measured by ric- The inner radius was 
fixed at 25 AU, the outer radius at 3 x 10'' AU, and the 
dust opacity was 0H5. The Xr i^) values for the 
SCUBA data show a clear minimum around ric = 1 x 10^ 
cm~^; the 1300 /im data favor ric = 1 — 3 x 10® cm^'^. 
The model with ric = 1 x 10^ cm~^ is a good compro- 
mise, and we adopt that value as our standard model in 
testing the effect of other parameters. We also explored 
the effects of using non-isothermal Bonnor-Ebert models, 
which may be a more realistic description of the density 
distribution. The differences in the predicted intensity dis- 
tributions were very small. 

The effects of chopping make the data quite insensitive 
to the outer radius. A model with an inner radius of 50 AU 
was essentially indistinguishable from the standard model, 
as was a model with Tq = 6 x 10^ AU. Decreasing the outer 



radius to 1.5 x 10'' AU has very httle effect on the SCUBA 
or 1300 /im fits. 

We have also considered the efi'ects of adopting different 
dust opacities. Models using 0H2 dust opacities differed 
little from those with 0H5 dust opacities in the fits to the 
SCUBA and the 1300 /im data, but the agreement with 
the SED worsened substantially. Because 0II2 dust has a 
higher opacity at long wavelengths, these models produced 
too much flux. The excess flux could not be decreased by 
varying Uc within the range favored by the intensity pro- 
files because the amount of mass in the beam changes little 
with ric- It is remarkable in fact that 0H5 dust opacity 
and a Bonnor-Ebert density distribution fits the SED so 
well. If we were sure that the cloud is hydrostatic and ther- 
mally pressure-balanced, so that a Bonnor-Ebert sphere is 
the correct physical model, then we could constrain the 
dust opacities quite strongly. 

The new far-infrared data (Ward-Thompson et al. 
2001) are about half the values previously reported, in- 
dicating lower Trf, and thus a lower value for the ISRF 
than assumed in our models. The last two models in Ta- 
ble 2 include these new data. The first retains the original 
ISRF, while the second scales the ISRF by a factor of 0.5, 
clearly improving the Xr for the SED. The latter model 
is shown in Fig. 5, including the temperature profile, the 
450/850 /im intensity ratio, the SED, and fits to the 450, 
850, and 1300 /im normalized intensity profiles. 

Given the substantial imccrtainties, the model fits the 
data well. In particular, the data suggest (with big uncer- 
tainties) a rise in the ratio toward the edge of the cloud. 
The model correctly predicts this behavior which arises al- 
most equally as a result of the increase in Tj, at the edge of 
the c;Ioud and as an effect of chopping. In principle, devi- 
ations of the data from the model could be used to study 
possible changes in dust properties with radius, but better 
signal-to-noise in the 450 /im data is needed before it is 
possible to make quantitative analyses. For the present, 
we note that any attempt to model the distribution of dust 
properties must incorporate careful considerations of spa- 
tial variations in the dust temperature and the efi'ects of 
chopping. We can also compare the actual values of the 
spectral index a to the observations (Shirley et al. 2000). 
The best model produces a = 2.1 in the 40" beam and 
a = 2.4 in the 120" beam (see Shirley et al. (2000) for 
the definition of a. These values are consistent within un- 
certainties with the a seen in the observations (2.0 ± 0.6). 

We may ask how our conclusions compare to the previ- 
ous studies of Andre, Ward-Thompson, & Motte (1996) 
(AWM). Firstly, AWM assume a distance of 160 pc and 
an opacity at 1300 /im of = 0.005 cm^ gm^^ of gas, 
whereas the corresponding value for the 0H5 dust opac- 
ity is Ki, = 0.009 cm^ gm~^ (assuming a gas to dust ratio 
of 100). In comparing results, we adjust the results of 
AWM to our assumed distance and opacity. The other 
main difference is that AWM assumed a constant = 18 
K, whereas we compute Td(r) and find it to be lower ev- 
erywhere, especially toward the center of the cores. Be- 
cause the lower Td depresses emission from the center sub- 
stantially (reducing the temperature from = 18 K to 
Td = 7.5 K results in a flux reduction by a factor of 4 at 
A = 1300/im), we find higher central densities (by factors 
of 3-10) and a smaller radius (by factors of 0.5-0.9) for 
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the flattened part of the density distribution. The small- 
est differences are found in the model with ric = 3 x 10^ 
cm~^. The mass inside 10^ AU is more similar; our mass 
(0.95 M0) is 1.4 times higher than that of Andre, Ward- 
Thompson, & Mottc (1996). Both studies agree that the 
density distribution outside the flat portion is well repre- 
sented by a power law (n(r) oc r~^) with p ~ 2. 

Bacmann et al. (2000) have re-analyzed the 1300 fim 
data, assuming a constant Td of 12.5 K and incorporating 
an analysis of the absorption at A ^ 7 fiia that the cloud 
shows in ISOCAM data. They also find a density dis- 
tribution that shows a distinction between a "flat" inner 
region and an envelope in which the density distribution 
can be described by a power law, roughly consistent with 
p = 2. They also find that Bonnor-Ebcrt spheres fit their 
data quite well. However, they find that the radius of the 
flat density core is 5000-6000 AU, 3 times larger than we 
find. This rfiat would correspond roughly to our Bonnor- 
Ebert sphere with ric = 1 x 10"'' cm^"^, which fits neither 
our SCUBA data nor the 1300 /im data, in our analysis. 

In principle, the absorption method does not depend on 
Td, thereby providing a complementary probe. However, 
the complication arises in the determination of the fore- 
ground emission, which veils the absorption by the core. 
Bacmann et al. (2000) use the 1300 /J,m data to constrain 
this foreground emission and conclude that the foreground 
is almost entirely zodiacal emission, so that the cloud must 
be illuminated entirely from behind. With this assump- 
tion, they find T(7/im) = 0.7; if there is more foreground 
veiling, the cloud could have r > 1 at 7 /xm, as our models 
indeed predict. Another variable is the opacity at 7 /xm. 
Bacmann et al. (2000) assume an opacity that translates 
to 3.1 cm^ gm~^ of gas, whereas the 0H5 dust opacity 
would imply 12.5 cm^ gm~^. Bacmann et al. (2000) also 
find evidence from the density profiles for an outer cut-off 
radius to the core. Since the core is flattened, this cut-off 
radius is variable, but taking the geometric mean of the 
two axes and correcting to a distance of 125 pc yields a 
radius of 2.6 x 10'' AU, similar to our outer radius, but 
outside the region that we effectively probe. While it is 
difficult to sort out all of these issues, it seems likely that 
the discrepancy between our models and those of Bacmann 
et al. (2000) imply that either they have underestimated 
the foreground veiling and hence T(7/im) or we have un- 
derestimated the dust temperatures. The new results in 
the far-infrared are well matched by our low Td, supporting 
our interpretation. Because the absorption measurements 
provide an independent measure of the distribution, this 
difference leaves a puzzle for further work. 

We conclude that Bonnor-Ebert spheres and 0H5 dust 
opacities reproduce the observations of L1689B very well, 
with ric ~ 1 X 10^ cm~"^. This density is higher and rfiat 
is smaller than have been found in other studies. 

7.2. LI544 

L1544 is at a distance of 140 pc (Ehas 1978). The 
SCUBA data have slightly better signal to noise than those 
of L1689B. While the source is clearly elongated (Shirley et 
al. 2000), we model it as a sphere. Since the radial profiles 
are spherical averages, this is the most (;onsistcnt compar- 
ison for the present. Ultimately, aspherical models should 
be constructed, but we address the likely effects in §8. In 



addition, we model 1300 /im data, subject to the same 
caveats as for L1689B: in particular, we do not expect to 
fit the intensity profile at large impact parameters. Mag- 
netic collapse models have been developed for this source 
(Ciolek & Basu 2000). Using their semi-analytic calcula- 
tion of dust temperature, Zucconi et al. (2001) find good 
agreement between this model and the data. Because we 
focus on spherically averaged data here, we do not attempt 
to compare this magnetic model to our data. 

The results are summarized in Table 3. First, con- 
sidering only simple Bonnor-Ebert spheres, a model with 
ric = 1 X 10^ cm~^ has the lowest Xr ^r the SCUBA data, 
as well as the 1300 /xm data and the SED. The model with 
ric = 3 X 10^ cm~^ has nearly as good a value of (S), but 
is slightly worse on the other measures. Density profiles 
appropriate to non-isothermal Bonnor Ebert spheres fit 
somewhat worse for ric = 1 x 10^ cm~^, but slightly better 
for Uc = 1 X 10^ cm~'^ than their isothermal counterparts. 

Most of the models did not include any far-infrared data 
in the calculation of the SED Xr value. The last two mod- 
els in Table 3 include the new far-infrared data in the SED 
Xr value. As with L1689B, the fit is greatly improved with 
a lower ISRF. The isothermal Bonnor-Ebert model with 
the ISRF decreased by 0.6 and ric = I x 10® cm~^ is shown 
in Figure 6. The spectral index, a, agrees well with the 
observations. 

Because models with very high Uc do not fit the data too 
badly, we have also considered the extreme limit of a sin- 
gular sphere {ric = 00), normalized to have the same den- 
sity at To as the Bonnor-Ebert sphere with ric = 1 x 10^ 
cm~^. This power law model has the density distribution 
of a singular isothermal sphere at a temperature of 10.4 
K (Shu 1977). We considered two models of this type: 
the first (iPL2) assumed a constant dust temperature of 
10 K. This model produced intensity profiles that are too 
strongly peaked in the inner core, confirming the conclu- 
sions of Ward-Thompson, Motte, & Andre (1999) that 
power law density distributions with constant Td cannot 
fit the data. A second model (tPL2) used the radiative 
transfer code to compute Td{r) self-consistently. Surpris- 
ingly, this model fits the SCUBA data reasonably well. 
While the fit to the 1300 data was clearly worse, our 
data cannot rule out a simple power law density model 
with a self- consistent Td{r) for L1544, given the current 
uncertainties in the data and the modeling. Figure 7 com- 
pares the results for the two power law models; it shows 
how important it is to include the self-consistent Td(r) in 
analyzing intensity distributions. 

As with L1689B, we compare our results to those ob- 
tained from the analysis of 1300 /zm emission (Ward- 
Thompson, Motte, & Andre 1999) and 7 /im absorption 
(Bacmann et al. 2000). The assumed distance is the 
same in all three studies, but (as for L1689B) the opaci- 
ties differ. Ward-Thompson, Motte, & Andre (1999) find 
He = 1 X 10® cm-3 [using the 0H5 value of ^,,(1300)], 
while Bacmann et al. (2000) find ric = 4 x 10^ or 1 x 10^ 
if 0H5 opacities are correct at 7 /xm. The radius of the flat- 
density-proflle inner core ranges from to 1600 AU in our 
models, compared to 2500 AU (Ward- Thompson, Motte, 
& Andre 1999) or 2900 AU (Bacmann et al. 2000). The 
latter authors also identify an outer (cut-off) radius of 8900 
AU. This is close to the size where our chopping and source 
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asphericity makes it difficult to constrain the density dis- 
tribution. As witli L1689B, our results arc qualitatively 
similar to other analyses, but differ quantitatively; the dif- 
ferences sTiggcst either that we have underestimated Td or 
that the optical depth at 7 /im has been underestimated. 

As for L1689B, Bonnor-Ebert spheres and 0H5 dust fit 
the data well. Our most striking conclusion that is differ- 
ent from the other studies is that our data do not rule out 
a singular power-law density distribution. 

7.3. L1512 

L1512 is also at a distance of 140 pc, but the submillime- 
ter emission is much weaker than from L1544. While the 
observational constraints are clearly much weaker than for 
the other sources, it is interesting to model L1512 because 
it seems to have a flatter distribution of intensity than the 
other sources. Shirley et al. (2000) presented only the 850 
/Ltm data, but the 450 /xm data does have a definable inten- 
sity distribution. Both the 450 /im and the 850 fjja data 
and some photometry constrain the models. The photom- 
etry was collected in tables in Shirley et al. (2000), except 
for a datum at 1300 yum from Ward- Thompson, Motte, & 
Andre (1999). 

The results for L1512 are given in Table 4. For this 
source we have only considered simple Bonnor-Ebert den- 
sity distributions — the poor data quality precludes the 
consideration of more sophisticated models. The best fit 
is for ric = 1 X 10^ cm~^, with ric = 3 x lO** cm""^ also 
acceptable. L1512 is clearly less centrally condensed than 
the other cores. 

The last two models in Table 4 include the new far- 
infrared data in the SED Xr value. The fit is clearly very 
bad for the full ISRF; the fit is much improved for an ISRF 
decreased by 0.3, but it is still worse than for L1689B and 
L1544. The last model is shown in Fig. 8. 

The problem with the SED is that the model fluxes are 
all higher than the observed values. However, good fits are 
obtained for the intensity profiles and the spectral index is 
consistent with the observations. Models with lower ric fit 
the lower fluxes better, but do not fit the intensity distri- 
bution. If the dust opacities are the same as in the other 
sources, fitting to the SED requires the core to be less 
massive. Bonnor-Ebert spheres calculated with Tk = 5 K 
contain less mass, and a model with ric = 3 x 10^ cm~^ flts 
both the SCUBA data and the SED best (Table 4). How- 
ever, such models are hard to justify on physical grounds; 
the radiative transport calculation produces Td{r) > 12 K 
everywhere (a model with Tk = 15 K fits the SED very 
poorly). While the density is too low to require Tk = Td, 
one would have to invoke unusually low gas heating to pro- 
duce Tk = 5 K. Another way to decrease the emission from 
the model cloud is to use different dust opacities; values 
between those of 0H5 and uncoagulated grains [e.g., col- 
umn 1 of Ossenkopf & Henning (1994)] would give about 
the right amount of emission. Such grains are plausible; if 
L1512 is a younger source, it will have spent less time at 
densities that are sufficiently high for grain coagulation to 
occur. 

L1512 was detected at 1300 /im by Ward-Thompson, 
Motte, & Andre (1999) but was too weak to map. It was 
also included in the absorption study by Bacmann et al. 
(2000), but the absorption was too weak to analyze. These 



results are consistent with our finding that the mass and 
central density are less than in L1689B and L1544. 

The radial profiles of L1512 can be matched if the source 
is a Bonnor-Ebert sphere with a central density smaller 
than those of L1689B and L1544. 

8. CAVEATS AND FUTURE WORK 

The fact that the masses of Bonnor-Ebert spheres, to- 
gether with the 0H5 dust opacities (Ossenkopf & Hen- 
ning 1994) match the submillimeter fiuxes so well suggests 
that both the physical and dust models may be reasonable. 
However, the opacity and mass effectively enter as a prod- 
uct; consequently, other mass and opacity combinations 
cannot be ruled out. 

Another caveat is that our models are one-dimensional. 
Some of the cores (L1544 most notably) are clearly as- 
pherical and (at least) two-dimensional radiative transfer 
models should be used to interpret the data. We can get 
a crude estimate of the effect of asphericity on the tem- 
perature distribution. Comparing T'd(r) for models with 
the same Uc but a factor of 2 smaller Vg simulates the ef- 
fects on heating of having one dimension smaller than the 
other. In fact, the differences are very small except near 
the surface of the cloud (< 0.2 K, or 3%, for r < 1 x 10'' 
AU). Semi-analytic calculations of in a two-dimensional 
geometry (Zucconi et al. 2001) find modest differences in 
the Td{r) along long and short axes. 

Differences in the assumed density law will produce dif- 
ferent Td{r) in general, as can be seen from Fig. 4. How- 
ever, if one plots the Td versus column density, the behav- 
ior is much more universal, as would be expected. This 
is shown in Fig. 9 for a variety of density profiles. The 
curves of Td{r) track each other until they approach their 
maximum value, where geometrical effects near the cen- 
ter of the cloud cause them to deviate from models with 
higher total column density. This figure may be useful 
for those modeling different density distributions who do 
not have access to a radiative transport code. For exam- 
ple, the T(j ~ 10 K computed from gray-body fits to the 
SED (Ward-Thompson et al. 2001) for L1689B and L1544 
would correspond to the model Td at about 1/3 the total 
column density into the cloud; in contrast, the radius at 
which Td = 10 K is about 4 x 10^ AU (Figs. 5 and 6), 
about 0.13 of the outer radius, enclosing 0.15 of the total 
mass. Clearly, gray-body fits to a constant Td weight ra- 
dius, column density, and mass differently, depending on 
the detailed density and temperature distributions. 

Local variations in the ISRF are another source of uncer- 
tainty. If local sources increase the strength of the ISRF, 
Td will increase. However, in practice, large changes in the 
ISRF strength arc needed for significant variations to oc- 
cur. Approximately, Td oc Li/(''+^), where L is the ISRF 
luminosity and [3 is the power-law index for the opacity 
oc v^^) (Doty & Leung 1994). For typical values of (3, 
the exponent is only 0.2 to 0.25, so that doubling the ISRF 
strength will increase the dust temperature by just 15- 
20%. Full radiative transfer models run with higher ISRF 
strengths confirm this simple analysis. Figure 3 shows the 
Td{r) calculated for the same physical model, but with the 
Black-Draine ISRF multiplied by a factor of 2 (except for 
the cosmic background radiation). The temperatures are 
higher at all radii by a factor around 1.14, as would be 
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predicted by the simple relations above. 

In fact, the new far-infrared data indicate a lower ISRF 
than we have assumed, even for L1689B, where the ISRF 
is generally considered to be enhanced. Thus lower Td 
are more likely than higher. Since Ward-Thompson ct al. 
(2001) describe their results in terms of enhanced ISRF, 
a word of explanation is in order. They use the ISRF of 
Mathis, Mezger, & Panagia (1983), which is considerably 
lower than the Black-Draine field that we use (cf. Fig. 
2a) in the mid-infrared to far-infrared region where they 
constrain the field. It seems that, for these cores at least, 
the correct value lies between the two standard choices 
for the ISRF. The far-infrared data strongly constrain the 
ISRF because they are exponentially sensitive to in 
this regime: future measurements and modeling will allow 
much better knowledge of the local ISRF. 

9. CONCLUSIONS 

We have modeled the emission from dust in pre- 
protostellar cores, including a self-consistent calculation 
of the temperature distribution (Td(r)) for an input den- 
sity distribution. We have also simulated the observations 
by convolving models with the observed beam and apply- 
ing chopping to the models. Using the calculated Td(r) 
has a substantial impact on the conclusions. Compared to 
earlier studies that assumed a constant Td{r), our models 
indicate smaller regions of relatively constant density (by 
factors of 0.5 to 0.9) and higher central densities (by fac- 
tors of 3 to 10). Indeed, for L1544, a singular, power-law 
density distribution cannot be ruled out. 

For the three sources we have modeled, there seems to be 
a sequence of increasing central condensation, from L1512 
to L1689B to L1544. While many more starless cores need 
to be modeled, it is possible that a new sequence of cores 
may be identified, in which increasing central condensation 
is the primary variable. The denser two of these cores, 
L1689B and LI 544, have spectroscopic evidence of con- 
traction motions (Tafalla et al. 1998; Gregersen & Evans 
2000), while L1512 does not (Gregersen & Evans 2000), 



consistent with this suggested sequence. 

It is interesting that Bonnor-Ebert spheres fit the data 
well, even though unstable spheres are needed in all cases. 
Magnetic fields may allow these nominally unstable ob- 
jects to persist, but in the absence of a suitable 2D/3D 
radiative transfer model we are not able to say whether 
or not the existing ambipolar diffusion controlled collapse 
models are consistent with the observations. 

Johnstone et al. (2000) have also used Bonnor-Ebert 
spheres to fit their SCUBA data in the p Ophiuchus molec- 
ular cloud. They find smaller outer radii and smaller de- 
grees of central condensation, but they have used models 
with constant Td = Tk- The smaller outer radii may re- 
flect the crowded conditions and higher ambient pressures 
in the p Ophiuchus region, compared to the "elbow room" 
available to the isolated cores in this study. More detailed 
analysis of the clustered regions, including calculations of 
Td{r), and more extensive studies of isolated regions will 
together delineate the differences in initial conditions for 
clustered and isolated star formation. 
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Fig. 1. — Plots of log density versus log radius in AU for the input physical models. In (a), Bonnor-Ebcrt spheres with constant kinetic 
temperature, Tk = 10 K and central densities from ric = 1 x 10'' to ric = 1 X lO'' cm~^ are shown. In (b), Bonnor-Ebert spheres with 
rtc = 1 X 10® cm^'^ arc shown for different, constant kinetic temperatures. In (c), a Bonnor-Ebert sphere with a variation in Tk, based on 
iteration with the radiative transport code, and assuming T/f (r) = T^[r), is compared to an isothermal {Tk = 10 K) Bonnor-Ebert sphere, 
with the same ric- In (d), a power law (PL2) corresponding to a singular isothermal sphere with Tk = 10-4 is compared to a Bonnor-Ebert 
sphere with ric = 1 x 10 cm~^. 
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Table 1 
Physical Models 



No. 


Type* 


IT— 




7^ — 7i — 


Tflat" 


IT- 

rio 


MS 






(AU) 


(cm^'') 


(K) 


(AU) 


(cni^'') 


(Mo) 


1 


BE 


3 X 10-* 


1 X lO" 


10 


16200 


1.84 X 10^ 


2.63 


2 


BE 


3 X 10^ 


3 X 10" 


10 


9300 


1.57 X 10^ 


3.15 


3 


BE 


3 X 10^ 


1 X 10^ 


10 


5110 


1.17 X 10^ 


3.06 


4 


BE 


3 X 10'' 


3 X 10^ 


10 


2990 


9.79 X 102 


2.75 


5 


BE 


3 X 10^ 


1 X 10*^ 


10 


1600 


9.38 X 10^ 


2.46 


6 


BE 


3 X 10" 


3 X 10*5 


10 


930 


9.90 X 10^ 


2.34 


7 


BE 


3 X 10^ 


1 X 10^ 


10 


500 


1.09 X 10''' 


2.34 


8 


BE 


6 X 10* 


3 X 10'" 


10 


2990 


2.34 X 10^ 


4.86 


9 


BE 


1.5 X 10'' 


1 X 10'^ 


10 


1600 


3.98 X 10^ 


1.40 


10 


BE 


6 X 10" 


1 X 10^ 


10 


1600 


2.50 X 10^ 


4.64 


11 


BE 


6 X 10" 


3 X 10*5 


10 


930 


2.74 X 10^ 


4.69 


12 


BE 


3 X 10" 


3 X 10" 


5 


6590 


6.59 X 10^ 


1.58 


13 


BE 


3 X 10" 


1 X 10^ 


5 


3600 


5.14 X 10^ 


1.43 


14 


BE 


3 X 10" 


1 X lO*' 


5 


1130 


4.82 X 10^ 


1.18 


15 


BE 


3 X 10" 


1 X 10^ 


15 


1970 


1.41 X 10^ 


3.82 


16 


tBE* 


3 X 10* 


3 X 10'" 


~7-15 


2500 


1.68 X 10^ 


3.55 


17 


tBE^ 


3 X 10" 


1 X 10^ 


~7-15 


1300 


1.68 X 10^ 


3.42 


18 


tBEi" 


6 X 10" 


1 X 10'5 


~7-15 


1300 


4.27 X 10^ 


7.25 


19 


tBE^ 


3 X 10" 


1 X 10^ 


~7-15 


400 


1.77 X 10^ 


3.44 


20 


iPL2 


3 X 10" 


1.56 X 10^ 


10.4 




1.08 X 10'' 


2.06 


21 


tPL2 


3 X 10" 


1.56 X 10^ 


10.4 




1.08 X 10^ 


2.06 



"BE is an isothermal Bonnor-Ebert sphere; tBE is a Bonnor-Ebert sphere with a temperature gradient; iPL2 is a power law, with an exponent 
of 2 and constant dust temperature of 10 K; tPL2 is the same power law with a dust temperature calculated by the radiative transport code. 
''The outer radius is To- 

'^Tic is the central density, at r = r;. 

"^Tn is the gas kinetic temperature used to compute the density distribution. 
"^Tfiat is the radius where the density drops to half the central density {ric)- 
f no is the density at the outer radius. 
®M is the mass enclosed within ro- 
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Table 2 
Models of L1689B 



No.'' 


Type 


ric 


M 




{5) 


Xr 


Xr 


Comments 






(cm^"^) 


(Mo) 


SCUBA 


SCUBA 


1300*^ 


SED 




3 


BE 


1 X 10-^ 


3.06 


16 


0.83 


1036-1560 


18 




4 


BE 


3 X 10^ 


2.75 


2.4 


0.38 


338-577 


17 




5 


BE 


1 X 10" 


2.46 


1.4 


0.26 


26-102 


16 




6 


BE 


3 X 10^ 


2.34 


5.1 


0.51 


30-22 


14 




7 


BE 


1 X 10^ 


2.34 


8.0 


0.73 


104-48 


11 




9 


BE 


1 X 10** 


1.40 


1.2 


0.26 


27-58 


6.6 


To = 1.5 X 10* 


8 


BE 


1 X 10^ 


4.86 


1.3 


0.26 


26-128 


18 


To = 6 X lO'^, n = 50 


5 


BE 


1 X 10*^ 


2.46 


1.4 


0.26 


25-102 


16 


ri = 50 


16 


tBE 


3 X lO'' 


3.55 


3.3 


0.38 


333-735 


20 




17 


tBE 


1 X 10^ 


3.42 


1.3 


0.26 


41 235 


19 




19 


tBE 


1 X 10^ 


3.44 


4.1 


0.52 


43 82 


16 




4 


BE 


3 X lO'^ 


2.75 


3.1 


0.45 


375 635 


212 


0H2 Dust 


5 


BE 


1 X 10^ 


2.46 


1.0 


0.27 


39-131 


237 


0H2 Dust 


5 


BE 


1 X 10** 


2.46 


1.4 


0.26 


26-102 


14c 


ISRF xl 


5 


BE 


1 X 10*^ 


2.46 


1.1 


0.27 


31-123 


1.5= 


ISRF xO.5 



"The model numbers refer to the physical model in Table 1. 

''The first value is for chopping by 120"; the second is for no chopping. 

°Only these two models include the new far-infrared data in the calculation of Xri the other models included previous far-infrared data. 



Table 3 
Models of L1544 



No.'' 


Type 


Tie 


M 


xi 


{5) 


xi 


xi 


Comments 






(cm^'^) 


(Mo) 


SCUBA 


SCUBA 


1300'= 


SED 




4 


BE 


3 X 10-^ 


2.75 


18 


0.63 


123 198 


2.1 




5 


BE 


1 X 10^ 


2.46 


0.73 


0.17 


14-26 


1.0 




6 


BE 


3 X 10*^ 


2.34 


3.0 


0.19 


118-49 


1.2 




7 


BE 


1 X 10^ 


2.34 


7.3 


0.39 


258-114 


2.0 




17 


tBE 


1 X 10'^ 


3.42 


2.8 


0.22 


28-88 


2.6 




19 


tBE 


1 X 10^ 


3.44 


2.0 


0.24 


160-88 


2.0 




20 


iPL2 




2.06 


22 


0.95 


480-247 


11 




21 


tPL2 




2.06 


3.0 


0.27 


149-210 


3.6 




5 


BE 


1 X 10** 


2.46 


0.73 


0.17 


14-26 


7.2^^ 


IRSF Xl 


5 


BE 


1 X 10'^ 


2.46 


1.1 


0.22 


13-31 


2.8^^ 


IRSF xO.6 



"The model numbers refer to the physical model in Table 1. 

''The first value is for chopping by 120"; the second is for no chopping. 

°Only these two models include the far-infrared data in the calculation of Xr- 
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Fig. 2. — In (a), the interstellar radiation field used in previous versions of the code ("MMP") is compared with that used here (Black- 
Draine). "MMP" is close to the ISRF used by Mathis, Mezger, & Panagia (1983), supplemented by a blackbody for the cosmic background 
radiation. The ISRF labeled Black-Draine uses the curve in Black (1994) for A > 0.36 /xm, and Draine (1978) for A < 0.36 /xm. In (b), the 
opacity per gram of gas (kv) for OH5 and OH2 dust, based on Ossenkopf & Henning (1994) and a gas to dust mass ratio of 100, is plotted 
versus wavelength. 
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Fig. 3. — Shows the eflfects on the dust temperature distribution of changing the ISRF from "MMP" to Black-Draine, and of using twice 
the average ISRF. A Bonnor-Ebert sphere with ric = 1 x 10® cm~^ and Tff = 10 K was used for all three calculations. 
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Fig. 4. — The density distributions for a series of Bonner-Ebert spheres are shown in (a). The resulting temperature distributions are 
shown in (b), with the same line coding; the densest models have the lowest central temperatures, in (c), the resulting intensity distributions 
at 450 iim arc shown. In (d), the results for 850 /im are shown. The same line coding is used for all panels, with the densest models showing 
the fastest fall-off with radius. The model emission in panels c and d has been convolved with the observed beam and chopping by 120" has 
been simulated, causing some of the drop around impact parameter, b ~ 10, 000 AU, for a distance of 125 pc. 
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Fig. 5. — A model for L1689B that fits the data well, a Bonnor-Ebert sphere with ric = 1 X 10 cm^"^ and a Black-Draine ISRF multiplied 
by 0.5. The temperature distribution is shown in (a). The observed SED is shown in (b) as circles with errorbars and an upper limit at 
90 /im; multiple values at the same wavelength arc data with different beams. The crosses are the model predictions for the same beams. 
The bottom two panels, (e) and (f), show for 450 and 850 /im, the observed, normalized intensity profile (circles and error bars) and the 
model (solid line), with simulated chopping. Panel (d) shows the data at 1300 fim of Andre, Ward-Thompson, & Motte (1996) and the model 
without simulated chopping (solid line) and with a simulation of chopping by 120" (heavy dashed line). Panel (c) shows the ratio of the 450 
and 850 fim normalized intensities with the same conventions. Note that normalization constrains the value to unity at the innermost point. 
At 125 pc, b = 10* AU corresponds to 80". 
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Fig. 6. — The same as Fig. 5, but for L1544. The model is a Bonnor-Ebert sphere with ric = 1 X 10^ cm and a Black-Draine ISRF 
multiplied by 0.6. The 1300 /im data is from Ward-Thompson, Motte, & Andre (1999). At 140 pc, b = 10"' AU corresponds to 71". 
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Fig. 7. — The normalized intensity distributions at 450 (a and c) and 850 (b and d) /im for a power law density distribution, n{r) = 
n{ro){r /ro)~^ , with n(ro) = 1.02 x 10'^, p = 2, and ro = 3 x lO** AU. The observations are of L1544. The top panels (a and b) show the 
model predictions with a constant T^{r) = 10 K, while the bottom panels (c and d) show predictions for Td(r) calculated self-consistently 
with the radiative transport code, using the full Black-Draine ISRF. 
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Fig. 8. — The same as Fig. 5, but for L1512. The model is a Bomior-Ebert sphere with ric = 1 X 10 ' cm and a Black-Diahie ISRF 
multiplied by 0.3. No data exist at 1300 /xm, so only the model predictions are shown (solid line for no chopping, dashed line for chopping by 
120"). At 140 pc, b = 10* AU corresponds to 71". 
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Fig. 9. — The dust temperature is plotted versus gas column density for a variety of models. The visual extinction is marked on the top 
axis, using standard conversions. The calculation was done for the full Black-Draine ISRF. 
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Table 4 
Models of L1512 



No.'' 


Type 


ric 


M 


Xr 


{5) 




Comments 






(cm-3) 


(Mo) 


SCUBA 


SCUBA 


SED 




1 


BE 


1 X 10" 


2.63 


3.2 


0.93 


6.6 




2 


BE 


3 X lO'^ 


3.15 


1.1 


0.62 


62 




3 


BE 


1 X 10^ 


3.06 


0.38 


0.41 


179 




4 


BE 


3 X 10^ 


2.75 


1.6 


0.59 


276 




12 


BE 


3 X 10* 


1.58 


0.31 


0.37 


5.1 


Tk=^ 


13 


BE 


1 X 10'^ 


1.43 


1.4 


0.57 


30 


Tk = 5 


14 


BE 


1 X 10^ 


1.18 


7.7 


1.2 


65 


Tk = 5 


15 


BE 


1 X 10^ 


3.82 


3.1 


0.73 


666 


Tk = 15 


3 


BE 


1 X 10=^ 


3.06 


0.38 


0.41 


ISS'^ 


ISRF xl 


3 


BE 


1 X 10^ 


3.06 


0.41 


0.43 


20'" 


ISRF xO.3 



"The model numbers refer to the physical model in Table 1. 

''Only these two models include the far-infrared data in the calculation of Xr- 



